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Abstract 

Wc show that a phcnomcnological hydrodynamic latticc-gas model of two-phase flow, devel- 
oped by Rothman and Keller in 1988 and used extensively for numerical simulations since then, 
can be derived from an underlying model of particle interactions. From this result, we elucidate 
the nature of the hydrodynamic limit of the Rothman-Keller model. 
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1 Introduction 

In 1986 it was discovered that certain mass and momentum conserving lattice-gas automata gave 
rise to the isotropic Navier-Stokes equations in the hydrodynamic limit In 1988, Rothman 

and Keller extended this discovery by introducing a hydrodynamic lattice-gas model of immiscible 
fluids Their model, and lattice Boltzmann variants thereof, have become an important tool 
for simulating the hydrodynamics of multiphase flow In both the original and RK lattice- 
gas models, the dynamics can be decomposed into two steps: In the first, the particles propagate 
along the lattice vectors to new sites; in the second, the particles entering each site collide by 
redistributing mass and momentum. 

In the Rothman-Keller (RK) model, the masses of the various immiscible fluid species and the 
total momentum are conserved locally, but the choice of collision outcome at each site of the RK 
model depends on the water-minus-oil order parameter, or "color," of the neighboring sites. See 
Fig. n for one set of possible collision outcomes that might be allowed by the conservation laws. 
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Figure 1: One Set of Possible Collision Outcomes for Immiscible Fluid LGA: The colors 
denote the particles' "color charge." The incoming state had one particle of each color and zero 
net momentum. Note that there are six possible outcomes that conserve these quantities. 

The flux of the color is determined for each such outgoing state, and its local gradient or field is 
determined by examining the neighboring sites. The negative of the dot product of this color flux 
and color field is then a measure of the propensity of outgoing particles to move to sites dominated 
by particles of their own type, and was called the color work by RK; their prescription was then to 
choose the outcome that minimizes this work in order to create cohesion and interfacial tension. (In 
case of a tie, the outcome is chosen randomly from among the states with minimal color work.) In 
later work, Chen, Chen, Doolen and Lee [H], and Chan and Liang noted that this minimization 
of color work is really just the low-temperature limit of a Boltzmann sampling procedure. 

This paper proposes a microscopic interpretation of the RK model. In the limit where the 
ratio of the mean-free path to the interaction range is small, we show that an arbitrary interaction 
potential can be expanded in terms of dot products of fluxes and fields, of which the RK model 
is merely the first term. This observation accounts for much of the success and utility of the RK 
model in describing the hydrodynamics of multiphase flow. 

2 Hydrodynamic Lattice-Gas Automata 

In lattice-gas models of hydrodynamics, incoming particles collide at each site x of a lattice C, 
in a manner to be discussed at length shortly, and then propagate to one of B neighboring sites 
x + Ci C, where i € {1, . . . , B}. (Note that some of the Cj's may be zero in order to accommodate 
"rest particles" in the model.) We suppose that the occupancy of each of these B channels can be 
represented by L bits (x) G {0, 1}, where i e {!,..., B} and i € {1, . . . , L}. Throughout the 
remainder of this paper, we shall illustrate various concepts by applying them to three concrete 
examples: 

• Example 1: In a lattice gas for a single-species Navier-Stokes fluid [SlIllIZli we take one bit 
(L = 1) in each direction that represents the presence or absence of a particle moving in that 



• Example 2: In a lattice gas for two immiscible fluids on the other hand, we might take 
two bits per direction (L = 2) so that there can be one bit for water particles n]^(x) and one 
bit for oil particles ?^P(x) in each direction. 

• Example 3: We consider the model of Chen, Chen, Doolen and Lee j^, in which there is 
one bit in each direction (L = 1), one "rest" particle direction i = R such that c/^ = 0, and 
only the rest particles feel an interaction potential. 

We suppose that there is a charge- like attribute gj(x) associated with the bits in direction i at 
site x, and we specialize to potential energies of the form 



direction. 
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where the factor of 1/2 prevents double counting. In Example 1, the charge-like attribute might be 
equal to the channel occupancy nj(x). In Example 2, on the other hand, the charge-like attribute 
might be the order parameter 

%(x)=nf(x)-nf(x) (2) 

which measures the excess of water over oil. In Example 3, we take 5i(x) = ni{x)5iR, since only 
the rest particles have a charge-like attribute. In what follows, we shall also have occasion to refer 
to the total (summed over directions) charge-like attribute of a site q{x) = J2iQi{^)- 

As noted above, any fluid model has some set of quantities that must be conserved in collisions. 
In Examples 1 and 3, we should demand conservation of the mass J2i^i{^) ^''^^ the momentum 

?T'j(x)cj. In Example 2, on the other hand, collisions must conserve water mass X^i^-l^C^)) oil 
mass J2i ^P(x) ^-iid total momentum X^J'^'i^lx) + nP(x)]cj. Alternatively, we can say that all three 
examples conserve total mass and the total momentum, while Example 2 also conserves the charge- 
like attribute, g(x). These conserved quantities naturally partition the set of all states of a given 
site(s) into equivalence classes of states with the same values for all of the conserved quantities. 
For example, Fig. ^ relevant to Example 2 above, illustrates the equivalence class comprised of the 
six possible collisional outcomes that may result when one water particle and one oil particle enter 
a single site on a two-dimensional triangular lattice from opposite directions (and, hence, with zero 
total momentum). 



3 Collisional Energetics 

We denote the postcoUision charge- like attribute with velocity Cj at site x by g-(x). Upon subsequent 
propagation, the charge ^^(x) will be at position x + Cj, and the charge q'j^x + y) will be at position 
x + y + Cj. This is illustrated in Fig. El The change in the potential energy due to both collision 
and propagation is then given by 

= EE + y)^ (ly + - C.I) - E E + y)<^ (lyl) 

= EE ^ii^m^ + y) (|y + c, - c,|) - (|y|)] 
x-y ij 

+ EE [(iii^Wji^ + y)-Qii^)Qji^ + y)] (^i\y\) 

x-y i,j 

= AVc + AVn, (3) 
where we have defined the contribution to AV due to the movement of the interacting particles, 

^ E E Qii^W^ + y) ['A (|y + c, - c,|) - <A (|y|)] , (4) 
x.y i,j 

and that due to nonconservation of the charge-like attribute, 

AVn ^ E Wi^Wi^ + y) - Qi^M^ + y)] <^ (lyl) • (s) 

Note that AVc vanishes for systems in which the interacting particles do not move, including our 
Example 3. Likewise, note that AVn vanishes for systems with a conserved charge-like attribute, 
including our Examples 1 and 2. 
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Figure 2: Change in Potential Energy: The change in the potential energy of interaction be- 
tween two charges moving in (possibly) different directions at (possibly) different sites is illustrated 
here. See Eq. ©. 

4 The Flux-Field Decomposition 

If the potential kernel (j){y) is analytic, we may expand AVc in a Taylor series in the ratio of the 
characteristic lattice spacing c to the characteristic interaction range y. This is done in Appendix^ 
where it is shown that every term of this series can be expressed as the complete inner product of 
a tensor flux with a tensor field. More specifically, we find that 

oo n r 

^^^ = EE E X{^)0£nA^), (6) 

X n=lr=\n/2'] 

where we have defined the (local) rth-rank tensor flux of outgoing particles, 

J-;(x)^^gKx)(^(g)c,^, (7) 

and where the (generally nonlocal) rth-rank tensor fields are defined in terms of (n — r)th-rank 
outgoing tensor fiuxes of neighboring sites, 

- (rS^(:)S'=»(y)OV^-,(x + y). (8) 

In all the above expressions, primes are used to denote dependence on postcollision values, 0^ 
denotes an r-fold outer product, O*" denotes an r-fold inner product, and we have defined the 
nth-rank completely symmetric kernel 

/2m— n \ /n—m ^ 



V I \ - 4>m{y) 

^n(y) = 2^ 7 rrper 

, , (n — my. 
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(9) 



^Since lattice gases are usually dense fluids for which particles undergo a collision at every step, c can also be 
thought of as a mean-free path. 
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where "per" indicates a summation over all distinct permutations of indices, and where we have 
defined the fohowing functions related to the derivatives of (p{y): 

ct>m{y)^[—r] m- (10) 



ydy 



The first few such kernels are thus given by 



[^i(y)]. = My)y^ (n) 
[^2{y)]i, = My)^ij + hiy)yiyj (12) 
[^3(y)]jjfc = My) (yi^jk + yjSik + yAj) + 4>3{y)yiyjyk (13) 

hiy) {yiyjSki + yiyk^ji + yiyi^jk + yjyk^a + yjyAk + ykyAj) + 
My)yiyjykyi- (i4) 

For completeness, we note that, since the zero-rank flux is just the charge, the portion of AV 
arising from nonconservation of the charge-like attribute may also be written in terms of these 
fluxes, 

= J E [^o(x) Jo'(x + y) - Jo(x) Jo(x + y)] My) (15) 

^ x,y 

as follows immediately from Eqs. © and 

Because the fields themselves depend on the post-collision fluxes at neighboring sites, it is 
generally necessary to include AVn in the collisional energetics. There are, however, some very useful 
exceptions to this rule. Some of the fluxes c7r(x) may be conserved, in which case J^i^) = c7r(x), 
since their values for incoming and outgoing states must be identical. If all of the fluxes that go 
into the calculation of a field £^ ,.{x) are conserved quantities, then we can write £^^^{k) = £n,ri^) 
as well. 

For example, let us return to consider our Example 2. At zeroth order, Eq. © indicates that 
AVn vanishes because is conserved. At first order (n = 1) the field £[ ^(x) depends only on the 
zero-order flux - namely the total order parameter (water minus oil) - at a given site, and this 
is conserved. Hence, the first-order energy change can be computed for the order-parameter flux 
for each outgoing state, and the field based on the incoming states of the neighbors x -|- y. If we 
restrict this set of neighbors to immediately adjacent sites, and look no further than lowest order 
in the Taylor expansion, we see that this reduces precisely to the RK model. The RK model is thus 
an approximation that is valid to first-order in c/y. While this fact may explain much of the utility 
of the RK model, in the following section we shall see that it may also indicate that the scaling 
limit of the RK model is more subtle than previously suspected. In any case, this also means that 
an alternative model that uses the exact potential energy, AV of Eq. Q, to sample post-collision 
states at each site will certainly be no worse than the RK model which is known to capture much 
of the phenomenology of immiscible fluid dynamics. 

Note that if both the flux and field in one of the terms of the Taylor expansion can be evaluated 
based on incoming quantities, then that term of AV is the same for all outgoing quantities, and 
therefore does not discriminate between them, so that it is necessary to go to higher order (at least 
until one encounters the first nonconserved fluxes) to obtain any interaction at all. We saw this 
with the vanishing of AT^i for systems with conserved charge-like attribute, such as Example 2 
above. To see this happen at higher order, let us consider the simpler Example 1 which, perhaps 
for this very reason, has been less studied. Both fluxes Jq and Ji are conserved, since they are 
the mass and momentum, respectively. It follows that the entire n = 1 term of AVc is the same 
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for all outgoing states. Nontrivial interaction between particles therefore does not even begin until 
the n = 2 term of the expansion, since J2 (which is, in fact, related to the pressure tensor) is not 
a conserved quantity. 



5 The Hydrodynamic Limit of the RK Model 

In order to obtain useful quantitative information from a hydrodynamic lattice gas, one must 
be careful to work in the correct asymptotic regime. This usually involves scaling the various 
dimensionless parameters of the problem with the Knudsen number Kn ^ A/L, where A is the 
mean-free path and L is the characteristic size. In incompressible Navier-Stokes flow, for example, 
one desires that the Mach number M = U/C, where U is the characteristic hydrodynamic velocity 
and C is the sound speed, scale with the Knudsen number M ~ ©(Kn). Since the viscosity z/ goes 
as the product of mean-free path and sound speed, AC, this implies that the Reynolds number Re 
scales as M/Kn ~ 0(ljl. We also demand that the Strouhal number St = Ut/L and the fractional 
density fluctuation 6p/pQ, where r is the mean-free time and po is the average background density, 
both scale as O(Kn^). This limit is well known [S] to reduce the compressible Navier-Stokes 
equations to their incompressible counterparts. Since, for a dense LGA, the mean- free path A 
goes as the grid size c, this means that in order to approach the continuum limit, every time one 
doubles the size of the lattice (halves Kn), one must quadruple the number of time steps (since St 
is quartered), and verify that the fractional density fluctuation (a measured "output" quantity in 
a lattice-gas simulation) is also quartered. Only when this scaling is verified can one be sure that 
one is working in the correct asymptotic regime. 

The presence of an interparticle potential adds an additional length scale ~ the range y of the 
force - and therefore a new dimensionless parameter X/y, or equivalently c/y. To derive the flux- 
field decomposition, we demanded that this ratio be small, but of order unity; that is, we did not 
scale this parameter with the Knudsen number. Operationally, this means that every time the size 
of the lattice is doubled, the range of the force in lattice units should be kept the same. If it is ten 
lattice units at one resolution, it should be ten lattice units at all resolution^. 

While Eq. (jlj for AVc is exact, the flux- field decomposition of Eq. © is usually used to ap- 
proximate AVc only to some specified order in c/y. Having determined that the RK model is 
just such an approximation, we are now in a position to examine how the error incurred by this 
approximation scales in the continuum limit. Let us compare the RK model to a variant of our 
model, in which we adopt the strategy of permitting the number of terms rimax that we retain in 
the Taylor expansion for AVc to increase in the hydrodynamic limit. We shall justify this strategy 
a posteriori. Specifically, let us take ^ more terms each time the lattice size N is doubled. If we 
take the system size in physical units L to be fixed, then the lattice spacing is c ~ L/N. For a 
dense lattice gas, the mean- free path is of order c, so the Knudsen number Kn scales as c/L. It 
follows that ^ 

n^s.^ = i log2 ^ =nQ-i log2 (Kn) , (16) 

where A^o ^ind no are constants. If we then take the interaction range in lattice units y/c to be 

^Note that this does not mean that the numerical value of the Reynolds number must be near unity. Rather it 
means that Re approaches a constant value in the scaling limit. 

''in some sense, this means that lattice artifacts never completely disappear, as the set of sites within a ten-lattice- 
unit radius of a given site are not distributed uniformly or isotropically. If this is deemed problematic, it may be 
possible to scale c/y with Kn~^^^, or in some other way such that both c/y and y/L vanish in the scaling limit; 
this has the attraction of completely removing such problems in the scaling limit, but further exploration of such 
considerations lies outside the scope of this paper. 
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fixed, then the error term e in the Taylor expansion goes as 

(^0"°Kn5i°S2(^). (17) 

For ^ = 1 (adding one more term to the series at each lattice refinement), it follows that by keeping 
the lattice spacing less than half of the characteristic interaction range, the error will scale as the 
Knudsen number; likewise, by making the lattice spacing less than a fourth of the characteristic 
interaction range, the error will scale as the square of the Knudsen number; and so on. One can also 
increase ^ to raise the power of the Knudsen number to which the error scales. Thus, for small hut 
finite values of c/y (order unity in the scaling limit), the error can be made to scale subdominantly 
to terms that are usually neglected in a Chapman-Enskog expansion anyway, providing the a 
posteriori justification promised above. 

A potential problem with the RK model is that it does not refine the definition of the energy 
in this way at each level of the scaling limit (.^ = 0), and so the corrections that it neglects may 
indeed matter in that limit. Of course, in most situations, one will choose a particular fixed value 
of nmax; in fact, all studies of the RK model to date have used rimax = l, simply because it has not 
been previously recognized that the RK model is only the first-order approximation to a more exact 
model. The RK model is known to exhibit certain anomalous phenomena; for example spurious 
currents are known to develop near interfaces, even if there is no bulk flow, and the surface tension 
is known to be slightly anisotropic [Oj. It is possible that such anomalies would be eliminated by the 
more exact treatment of the scaling limit advocated here, but a numerical test of this conjecture is 
outside the scope of the present work. 

Of course, such anomalies may be regarded as tolerable as long as one appreciates that one 
is working only to lowest order of an asymptotic series. Indeed, because the series is asymptotic, 
there is no point in being overzealous about the value of nmax? since the series may begin to 
diverge at some point. There is at least one situation, however, in which the observation that 
it is necessary to let nmax scale with Knudsen number may be critically important, and that is 
when one is studying the scaling of (possibly divergent) quantities with system size. For example, 
Rothman and Flekk0y JHI recently studied the scaling properties of fluctuating interfaces using the 
RK model, measuring, among other things, the saturated width of the interface as a function of 
system size. Superimposed upon the usual power-law behavior of the saturated width, they found a 
logarithmic correction that resisted theoretical explanation. We suggest that this anomaly may be 
due to scaling the system size with fixed Jimax (since the RK model effectively fixes nmax at unity), 
rather than letting nmax increase linearly with system size as advocated here. Again, it would be 
interesting to verify this conjecture by redoing the numerical experiment using our model, but that 
is outside the scope of the present work. 

The main point of this section has been to demonstrate that the flux-field decomposition allows 
one to understand in what sense the RK model is an approximation to an exact interaction potential, 
why it may be used for trial state sampling, how it might be corrected at higher order in c/y, and 
why its scaling limit may be more subtle than previously suspected. We emphasize that the strategy 
of scaling n^nax with Kn was invoked only to facilitate this discussion of the scaling limit, and we 
are certainly not suggesting that it be used in practical simulations. If one is planning to work 
at any order of c/y at which AV^ involves nonlocal interactions, it makes much more sense to use 
Eq. ((3} directly than to try to deal with higher order terms of the fiux-field decomposition. For 
Monte Carlo sampling of the outgoing states, it is possible that the lowest-order local term of the 
flux-field decomposition could be used for sampling, while the exact expression for AT4 could be 
used for the acceptance criterion. 
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6 Conclusions 



We have shown that the lattice gas model developed by Rothman and Keller for immiscible-fluid 
hydrodynamics can be derived from an underlying model of particle interactions. From the en- 
hanced understanding provided by our observation, we elucidated the nature of the hydrodynamic 
limit of the Rothman-Keller model, demonstrating that it is more subtle than previously suspected. 
Though practical simulations of the particulate model are likely to be significantly more compute- 
intensive than the original version of the Rothman-Keller model, this work is off'erred in the spirit 
that it is always useful to know the exact model corresponding to any given approximation. We 
hope that this work helps to provide some theoretical basis for these models' success, and perhaps 
for their ultimate improvement. 
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A Derivation of Flux-Field Decomposition 

To derive Eq. for AVc, we begin with the Taylor expansion of a function of the magnitude of a 
displaced vector, 

in I |\ ^ n.cpniyy) I \lm.—n I \n—m / ^ o\ 

^(|y + -l) = E r! E o»-^(2^-n)!(n-m)! (y-"^ ' ' ^^^^ 

n=0 m,= \n/2\ ^ ' ^ ' 

where e has been introduced as an expansion parameter to keep track of the order in c/y (it is 
numerically equal to one), y = |y|, and we have defined the following functions related to the 
derivatives of 4>{y): 

(19) 

We let c — > Cj — Cj and use the binomial theorem to write 

/ \2m—n I \2m—n 

2m~n /r) \i 



and 



(c-cf-™ = (|c,f + |c,|2-2c,.c,)' 



n-mn-m-k fr. _ Ic • P'^ Ic • 

f^, k\p\{n - m - k - p)r ' ■ ^ > 



Inserting these into Eq. @ for AVci we get 

Ay, = 

^EE'^Kx).Kx + y)E^ E 

x,y i,j n=l ■ m=[n/2] 

2m-n n-mn-m-k / -, n „_m+i-fc-p^ | 
\ ^ \ ^ \ ^ \ — ^) 

f-i f-^r. n 2^+n\k\p\{2m - n - ly.in - m - k - p)\ 
l=Q fc=0 p=0 ^ ' ^ ' 

<Pm{y) (y • c,)'""""' (y • c,)' \c,nc,\'^ (c, • c,)"""-'-^ . (22) 

Eliminating p in favor of the new summation index r = n — m + l + k — p, and reordering the 
summations, this becomes 

AK = 

oo n n n 2m-n min(n-m, [{r-/)/2j ) 

iE^EEE^^(xkKx + y) EE E 

n=l r=0 x,y ij' m=[n/2] /=0 fe=max(0, m-n+i — /) 

.\2k\ ,\2n-2r-2m+2l+2k 



-lYn\<t>miy)\cir\c 



2n-r-m+i+2ki\j^\(^n - r - m + / + k)\{2m - n - /)!(r - / - 2k)\ 
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where we have adopted the convention that a sum is zero if its upper hmit is less than its lower 
limit. By reinterpreting dot products raised to the s power as the s-fold inner product of two s-fold 
outer products, we can rewrite this as follows 

n=l ■ r=0 \ / x,y ij 



n—r /n—r 



(24) 



where 0^ denotes an r-fold outer product and denotes an r-fold inner product, and where we 
have defined the kernel 

n 2m-n mm{n-m,[{r-l)/2]) 

^n(y) -EE E 

m=\n/2~\ 1=0 fc=max(0,m— n+r'— 

r!(n - r)l(pm.iy) 



2n-r-m+i+2fc;!/,|(-^ - r - m + / + ky.{2m - n - l)\{r - I - 2k)\ 

/2m— n \ /n~m ^ 



y — 

on- 



y j ^ 09 1 

n\(l)m{y) 



m=\n/2'\ 



''(2m — n)!(n — m)! 



/2m— n 



(25) 



In the very last step above we performed the sums over k and /. 

The expression for the potential energy, Eq. H24|) . has a remarkable symmetry, inherited from 
Eq. (121). By making the substitutions 



I 
j 

X 

y 



n — r 

3 
i 

X + y 

-y 



(26) 



and noting that /Cri(— y) = (— l)"'^n(y)) we can see that the rth term of Eq. H24() is equal to the 
(n — r)th term. It also follows that the kernel Kn can be chosen to be completely symmetric under 
interchange of any two of its n indices. If we notice that the combinatorial factor 



(27) 

2"-™(2m - n)! ^ ' 

is precisely equal to the number of distinct ways to assign n indices to the tensor ((^^"^~" y) (g) 
({^""'"l), and recalling that only the symmetric part of /C.„ matters, we can rewrite the kernel, 
Eq. (|25|). in the remarkably compact form of Eq. Q. The simplicity of this result suggests that 
there may be an easier way to derive it. 

If we now introduce the completely symmetric rth-rank outgoing tensor fluxes, as defined in 
Eq. (I?!), then Eq. pi]) may be written as 



oo n "f^ / \ ^ n—r 

i E E 3 E(-ir ! Xi-) O ^n(y) O ^Lri- + y)- 



2 ^ ^ n! 

x,y n=l r=0 



(28) 
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Because of the symmetry, Eq. ()26|). we can remove the factor of 1/2 and sum r from [n/2] to n, 
instead of from to n. The exception to this arises when n is even and r = n/2; in that case, 
the factor of 1/2 must be retained. We accommodate this case by dividing by 1 + 6n,2r, where 5 is 
the Kronecker delta. If we then define the rth-rank tensor fields, as in Eq. (jHJ, the expression of 
Eq. Q for the potential energy change follows immediately. 

As noted in the text, the primes on the fluxes and fields indicate that these are evaluated using 
the outgoing (post-collision) states. Thus, at each order n, Eq. @ expresses the change in potential 
energy due to the propagation step as the sum of n terms, the rth of which is an r-fold inner product 
of an rth-rank outgoing tensor fiux, Eq. (jT)), with an rth-rank tensor field, Eq. 

Note that we can rearrange the order of summation of n and r to write Eq. © in the following 
alternative format 



This makes it clear that at each order in r a new rth-rank tensor flux must be introduced, and the 
corresponding rth-rank fleld is the sum of r terms. The disadvantage of writing the result in this 
way is that each term of the outer sum over r contains terms of differing order in e. 

Finally, we note in passing that tensor fluxes and fields have been considered in the context 
of a LGA model of microemulsions While that model did indeed derive its update rule from 
considerations of single-particle interactions, it did not employ the method used here. In particular, 
where Cj — Cj appears in Eq. (jH)), that study took only — Cj, and the expansion was not carried out 
to all orders. On the other hand, that work also included vector-valued charge attributes, in 
order to properly model the orientation of the surfactant molecules. It is likely that a more exact 
formulation, analogous to our Eq. (jSJ, exists for this microemulsion model, and may well connect 
this work with the static lattice models of microemulsions due to Matsen and Sullivan ^2] • 




(29) 
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